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Using a recently introduced algorithm for simulating percolation in microcanonical (fixed- 
occupancy) samples, we study the convergence with increasing system size of a number of estimates 
for the percolation threshold for an open system with a square boundary, specifically for site per- 
colation on a square lattice. We show that the convergence of the so-called "average-probability" 
estimate is described by a non-trivial correction-to-scaling exponent as predicted previously, and 
measure the value of this exponent to be 0.90 ± 0.02. For the "median" and "cell-to-cell" estimates 
of the percolation threshold we verify that convergence does not depend on this exponent, having 
instead a slightly faster convergence with a trivial analytic leading exponent. 



I. INTRODUCTION 

Percolation Q is one of the most fundamental and 
widely studied systems in statistical physics. Theoretical 
studies of percolation models and applications of percola- 
tion theory to physical systems have spawned thousands 
of papers over the last few decades. Even so, there are 
some substantial gaps in our understanding of percola- 
tion. For example, we have at present no exact value for 
the position p c of the percolation threshold for site per- 
colation on that simplest of two-dimensional lattices, the 
square lattice. And in three dimensions we have almost 
no exact results whatsoever. Because of this, numerical 
methods have played an important role in the study of 
percolation. In this paper we consider a class of meth- 
ods for estimating p c for site percolation using finite-size 
scaling, and show how various estimates of p c in this class 
scale with varying system size in two dimensions. 

The methods studied here for measuring p c are widely 
used and are all based upon consideration of the cross- 
ing probability function Rl(p) El [§■ This function 
gives the probability that a connected path crosses the 
system from one boundary segment to another, at site oc- 
cupation probability p and system size or length-scale L. 
Some examples of these estimates are: 

1. The renormalization- group (RG) fixed-point esti- 
mate pkg{L) = P*(L), where p* is the solution to 
the equation || [| 

R L ( P )=P (1) 

2. The average value of p at which crossing first oc- 
curs § §: 

p av (L) = ( P ) = f P R' L (p)dp 
Jo 

= 1- f R L (p)d P , (2) 
Jo 



where the last equality follows from integrating by 
parts. The prime indicates differentiation with re- 
spect to p. 

3. The estimate p/j e (L) corresponding to the point 
where Rl(p) equals its universal infinite-system 
value R c = Rqc(Pc) (which is determined by the 
system shape and boundary conditions) [|| : 

R L (p) = Rc- (3) 

For a square system, where R c = |, this estimate 
Po.b{L) corresponds to the median of the distribu- 
tion R'{p). A related estimate ph+ v {L) for rectan- 
gular systems is the value of p at which the hor- 
izontal and vertical crossing probabilities sum to 
unity § 0, | : 

R ( L h) (p) + R^(p) = l. (4) 

This estimate is identical to pn a (L) — po.s(^) when 
the boundary is a perfect square. 

4. The estimate Pmnx(L) which is the value of p at 
which R' L (p) reaches a maximum (or equivalently, 
where Rl(p) is at its inflection point) 

R'l(p) = 0. (5) 

5. The cell-to-cell RG estimate, which is the point 
where two systems of different size have the same 
value of R g . One possible choice for this estimate, 
p^} c {L), is the value of p at which 

R L {p) = R L -i(p), (6) 

while a second, p^_ c (L), is the point at which 

R L (p) = R L /2(p)- (?) 
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In order to use these estimates to determine the thresh- 
old precisely, wc need to know the manner in which they 
converge to p c as L — > oo. While it is possible to sim- 
ulate very large systems for which finite-size effects may 
be quite small, the statistics for such simulations are still 
relatively poor because of the small number of samples 
that can typically be generated. In most cases, better re- 
sults can be derived by doing more simulations on smaller 
systems, and this requires that the finite-size behavior is 
characterized accurately. 

It is usually assumed (based upon very general scaling 
arguments) that all finite-system estimates of the perco- 
lation threshold converge to the bulk value p c as 

PvM-p^cL- 1 '" (8) 

where c is a system-dependent constant, and v is the 
exponent governing the correlation length £, such that 
£ ~ \p ~ Pc\ ~ v ■ (For the two-dimensional systems we 
will be looking at in this paper v = |.) Well known 
exceptions to this behavior are a few highly symmetric, 
self-dual systems, such as bond percolation on a square 
lattice with a square boundary, and site percolation on a 
triangular lattice with a rhomboidal boundary; in both 
these cases Rl(p) is perfectly symmetric about P = \ for 
all L and all the estimates above give p c = \ exactly. For 
these systems, one could consider the constant c above 
to be zero. 

In Ref. however, it was argued that for non-self-dual 
systems with a square boundary, such as site percolation 
on a square lattice (where because of the non-duality 
the estimates show finite-size effects) , the convergence of 
most of those estimates is faster than given by Eq. (|^). 
This is an observation of some practical significance, since 
this particular system (site percolation on a square lattice 
with a square external boundary) is one of the most com- 
monly studied systems in percolation. Similar arguments 
also apply to other symmetric two-dimensional crossing 
problems, such as a system with a rhomboidal boundary, 
which is commonly used when simulating triangular and 
honeycomb lattices. 

The arguments of Ref. ^| were based upon the hypoth- 
esis that 

Rl(p)~ h(x)+L- 1 f 1 (x) + ... (9) 

for large L, where x = (p — p c )L 1/>u , fo(x) represents the 
universal part of R, and fi(x) represents the first-order 
correction to the scaling limit. The choice of L^ 1 as the 
leading order of the correction was based on numerical 
measurements of R at p c , and can be derived from the 
assumption that the system is effectively slightly rectan- 
gular in shape, because of the different types of boundary 
conditions applied along the two principal axes |J. For 
small x, it was assumed that 

fo(x) — a + aix + a 3 x 3 + ... (10) 
fi(x) = b Q + b lX + b 2 x 2 + ... (11) 



where clq = \ by the symmetry and self-duality of the 
square boundary. The same symmetry also implies that 
fo — ^ is an odd function in x and hence that a n = for 
all even n > 0, as above. (To see that f — \ is odd, note 
that for the perfectly dual system of bond percolation 
on a square lattice, R — | is an odd function of x for 
all L, including L = oo, and by universality systems with 
other underlying lattices must behave the same way in 
the scaling limit.) In Ref. ^, no particular assumption 
was made about the behavior of fi(x), other than its 
analyticity about x = 0. 

Note that the form of fo(x) in Eq. (|To| ) is not entirely 
universal, because the independent variable x should in- 
corporate a metric factor which depends upon the under- 
lying lattice @ 0. For convenience, however, since we 
are considering only the one system of site percolation 
on a square lattice in this paper, we do not include that 
factor here. 

Once the above assumptions, Eqs. (|l(]) and ( pT| ) are 
made, the convergence of the various estimates of p c is 
straightforward to analyze, and one finds that while the 
RG estimate does indeed converge according to Eq. (||) (a 
result that has been verified in many numerical studies), 
the rest of the estimates above should converge according 
to the faster behavior 

Pest {L)-p c ~cL- x - x l\ (12) 

where the constant c varies from estimate to estimate. 
Simulations reported in Ref. for systems of size up to 
1024 x 1024 sites verified this convergence for the esti- 
mate po.5 to high accuracy. The estimates p c -c, Pmax 
and p av were studied in Ref. |5| using only exact enumer- 
ation results for systems of sizes up to 7 x 7 which give 
polynomials for R (see the appendix), and while the be- 
havior of these results was found to be roughly consistent 
with ( |l2"| ) , the uncertainly due to higher-order corrections 
was large. 

Following the publication of Ref. [|, Aharony and 
Hovi |h], [l^] argued that the irrelevant scaling variables 
in the renormalization-group treatment of percolation 
imply a slower leading-order convergence of Rl to its 
infinite-system value, characterized by an exponent to, 
whose value was deduced from the Monte-Carlo work of 
Stauffer p]| to be about u> — 0.85. (Note that Aharony 
and Hovi used Q\ to denote the exponent we call ui.) A 
variety of series expansion results from the early 1980s 
were also analyzed to give values for this exponent rang- 
ing from 0.89 to over 1 Q, ||. 

The argument given by Aharony and Hovi implies that 
the leading terms in the expansion of Rl{p) should in fact 
be 

Rl{p) ~ /o(a:) + L-"f w (x) + L- l h{x) + ... (13) 
where 

f u (x) = ax + c 2 x 2 + . . . (14) 
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Hovi and Aharony argued that the constant term cq is 
zero for a square system, because at p c there are no cor- 
rection terms for the square-lattice, bond-percolation sys- 
tem, and corrections due to irrelevant variables should be 
universal. They also argued that fi(x) should be even, 
so that b\ = in ( |il"l) , also by symmetry and self-duality. 
They discussed various consequences of these assump- 
tions, and presented numerical evidence that the term 
containing the exponent u> is indeed the leading correc- 
tion term, by showing that the behavior of Rl{p) for 
large x (that is, for p ^ p c ) was better fit with such a 
term than without it. However, the procedure they used 
did not allow them to determine the value of u accurately, 
because of the numerical difficulty of finding Rl (p) pre- 
cisely for all p. Furthermore, they did not study the 
convergence of the other estimates given above. 

Recently (if], |l7j , the present authors have shown that 
quantities like Rl(p) can be studied efficiently for all p 
by first finding the crossing probability Rl.u in a "micro- 
canonical" system of exactly n occupied sites, and then 
convolving with a binomial distribution to derive results 
for the corresponding "canonical" system thus: 

R L (P)=J2( ) P n (l-p) N - n RL,n, (15) 

where N — L 2 for site percolation on a square lat- 
tice. The microcanonical crossing probability is found 
using an efficient cluster-joining algorithm employing 
data structures based on trees and a fast method is em- 
ployed for checking for percolation on the fly during the 



progress of the calculation. (While many of the ideas 
incorporated in this method were put forward previ- 
ously @ J^, ||, ||, this was the first time that 
all of these components were combined in this way, for 
the purpose of finding R efficiently.) In Rcf. [l6| we stud- 
ied the function corresponding to R for the probability 
of a cluster wrapping around the boundary of a periodic 
system on a torus. In the present paper, we describe how 
that method can be implemented for crossing of an open 
system, and we report results from some large-scale sim- 
ulations. The results allow us to determine accurately 
the behavior of all of the estimates above, and to test 
the theoretical predictions that follow from Eq. (|l3"|). As 
we will see, the appearance of the "irrelevant" term in 
the scaling of some estimates is confirmed, and a new, 
more precise value of u> is found. 

The outline of the paper is as follows. In Section || 
we derive the expected scaling behavior of the various 
estimates of p c , assuming the form (|Hf). In Section |l| 
we describe our numerical method, and in Section |M we 
present the results of our simulations. In Section we 
give our conclusions. 

II. CONVERGENCE OF ESTIMATES 

If we assume Eq. (|l^) to be a correct description of 
the behavior of Rl(j>), it is straightforward to deduce the 
resulting convergence of the various estimates for p c . In 
the following paragraphs we derive the leading correction 
term for each estimate, or leading two terms when their 
powers are close to each other. 



I 

1. The RG fixed-point: The relevant terms of Eq. (nh are 



which implies 



Pkg{L) =p c + 



i + a x x + a 3 x 3 + . . 

a 3 (p c - 



xL- x l\ 



Pc 



1)3 



a i 



+ 



L" 1 / 1 ' + 0(1/ 



■2/vs 



(16) 
(17) 



The term in brackets is the value of x that is the solution to fo(x) = p c . 
2. Average-probability estimate: Eq. (^) gives 

p av (L) = 1 - L~ x l v f 1 [f (x) + L-«f u {x) + L- l h{x)] dx, 



(18) 

where xq — ~p c L 1 l v and x\ = (1 —p c )L 1 l v , which are the values of x at p = and 1 respectively. Noting that, 
since fo(x) is odd about /o(0) = \ and approaches 1 as x — > oo, we have 



L~ x l v \ fo(x)dx^L- 1 ^x 1 = l-p c 



and we then get 



Pa 



{L)=p c + L-"- l t v f^dx + L- 1 - 1 ^ fi(x)dx, 



(19) 



(20) 



where we have extended the limits of the integrals to ±oo. This result is also implied by Eq. (40) of Ref. [lj, 
for n = 1. The order of the next correction depends upon the higher-order corrections to (([3]). 
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3. Median-p estimate: Eq. (||) gives 

| + a 1 x + b Q L~ 1 + Cl xL~".. . = ±, (21) 

which implies 

PoAL) =Pc- —L~ 1 ~ 1 / 1 ' + CKZ- 1 — V"). (22) 
ai 

4. Maximum estimate: Eq. (||) gives 

6a 3 x + 2b 2 L~ 1 + (2c 2 + 6c 3 x)L- u + . . . = 0, (23) 

which implies 



/ r \ c 2 t-oj-1/v ^2 J ._i_i/ J/ ri / r -2o;-l/i/s 
Pmax(L) = p c - L 1 — - L ' +U{L * ) 

3a 3 3a 3 



5. Cell-to-cell estimate: Eq. (H) gives 



which implies 



(2) 

Likewise, for p c _ c we have 



(2) 
Pc-c 



(24) 



a + a x (p - p c )L 1/,y + feoi" 1 + ci(p - p e )lM"- u + . . . 

= a + ai (p- Pc )(L - \fl" + b Q (L - I)" 1 + d(p - p c )(i - I) 1 / 8 — + . . . (25) 



p£> (L) = p c + — zaL- 1 - 1 ^ + □(L- 1 -^- 1 ^). (26) 
ai 



(£) - + ai(1 _ & ° 2 _ 1/ , ) ^ 1 ' 1/t/ + Ofi- 1 -- 1 /"). (27) 



Thus, the scaling of prg, Po.5, an d p c -c is unaffected to 
leading order by the presence of the exponent u>. How- 
ever, p av and p max are affected by u> to leading order, scal- 
ing as L~ u ~ x / V , slightly slower than predicted in Ref. ^ 
(order L~ 1 ~ 1 / u ). Note that b\ does not enter in any of 
these results, so the leading scaling does not change if b\ 
is equal to zero, as argued to be the case in Ref. O. 



III. PROCEDURE 

We have performed simulations to test the scaling hy- 
potheses above using the algorithm described in Refs. [l6| 
and Briefly, sites are occupied one by one in random 
order starting with an empty lattice. Occupied sites form 
contiguous clusters, each of which is identified uniquely 
by the site label of a chosen single site within the clus- 
ter, which we call the "root site." Other (non-root) sites 
within a cluster possess pointers that point either directly 
to the root site, or to other sites within the cluster such 
that by following a succession of pointers one can get 
from any site to the root. A newly added site is consid- 
ered to be a cluster of size one, which is its own root site, 



and bonds are then added between it and any adjacent 
occupied sites. The clusters to which sites at either end 
of such a bond belong are identified by following pointers 
from them to their corresponding root sites, and if the 
root sites found are different we conclude that two differ- 
ent clusters have been joined by the addition of the bond. 
We represent this by adding a pointer from the root site 
of one of the clusters to the root site of the other. Smaller 
clusters are always made subclusters of larger ones, and 
all pointers followed are subsequently changed to point 
directly to the root of their own cluster. The net result 
is an algorithm that can calculate Rl,u (and many other 
observable quantities) for all values of n in average run- 
ning time which is of the order the area of the lattice, or 
0(L 2 ) for a square lattice. 

In our previous calculations using this algorithm we 
measured the probability of the existence of a cluster 
that wraps around the periodic boundary conditions of a 
toroidal lattice. In this paper we are interested instead in 
the existence (or not) of a cluster that spans an open sys- 
tem along one given direction. There are (at least) three 
efficient methods for detecting spanning of this kind, two 
of which are described in detail in Ref. [L7l and all of which 
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we have used in the present work. In the first method, we 
use the same pointer-based trick that we used in Rcf. [l6| 
to detect wrapping with periodic boundary conditions, 
but start out with an (L+ 1) x (L + 1) lattice in which one 
horizontal row of sites is fixed to be permanently empty 
and one vertical one is fixed occupied. Occurrence of a 
wrapping cluster in such a system is then exactly equiva- 
lent to occurrence of a spanning cluster in the horizontal 
direction in an open system with dimensions L x L. 

In the second method, two complete rows of sites at 
the top and bottom of an open (L + 2) x (L + 2) lat- 
tice are fixed permanently empty, and two columns of 
L sites on the left and the right sides of the lattice are 
fixed occupied. The two columns of occupied sites form 
two initial clusters on the lattice. By following pointers 
from one site in each of these clusters it is then simple 
to determine whether the two clusters have been amalga- 
mated by other sites added between them: if they have 
the same root site they have been amalgamated, oth- 
erwise they have not. Performing this check after the 
addition of each site to the lattice, we can detect when a 
spanning cluster on the L x L open lattice first appears. 
The two methods have comparable running times and 
give compatible results. The second is somewhat simpler 
to implement. 

In the third method, which was used for the major- 
ity of the simulations here, we consider an L x L open 
lattice and keep track of the minimum and maximum 
x- and y-coordinates for sites in each cluster, updating 
their values as necessary when clusters are joined. When 



percolation threshold are then evaluated directly accord- 
ing to Eqs. (||), (||), and (^). The remaining estimate, 
p av , could be found directly by performing a numerical 
integral over Rl(p), but a better method is to use the 
following exact formula: 

Pav = 1 - / Rh(p) dp 

Jo 

= l-p^JRL,nJ o 1 p n (l-p) N - n dp 
1 N 



(29) 



n=0 



Using Eq. ( |28|) this can also be written directly in terms 
of P L ,« as 

N n N 
n=Q n'—O n=0 

(30) 

which means that the canonical average position of the 
percolation threshold is N/(N + 1) times the microcanon- 
ical average (1/N) J2nPL, n and no convolution is neces- 
sary to find its value. Higher moments of the distribution 
R' can be found in a similar fashion. For the second mo- 
ment, for example, we have 



ecessary when clusters are joined. When pi 
L — 1 for a cluster, we know that the (p 2 ) = / p 2 R' L (p) dp 

p lattipp in thp hnri'znntnl Hirpr'tinn nnH "'0 



cluster spans the lattice in the horizontal direction, and 
similarly for vertical crossing. This method allows one 
to check for both crossing events simultaneously, and it 
is also efficient and easy to program. A similar method 
was used in Ref. [24| for simulations of the Ising model. 

Since in the present calculation we are only interested 
in the existence or not of a system-spanning cluster, we 
can stop the simulation once a spanning cluster is de- 
tected, as spanning must also occur for all higher values 
of n. This produces about a 40% saving in running time. 
Each simulation then produces just a single number, the 
value of n at which a spanning cluster first appears (or 
two numbers if we check for spanning in both the hor- 
izontal and vertical directions). Making a histogram of 
these values over many runs of the algorithm, we derive 
an estimate of the probability Pl,u that the system first 
percolates when the number of occupied sites reaches n. 
This probability is related to the desired function Rl^ 
according to Pl.u — Rl,u — Rh.n-i, and hence 



N 



(N + l){N + 2) 



1 



JY 



J2n(n + l)P L , n . (31) 



(N + l)(N + 2)^ Q 

To find where R'l(p) — for the estimate p max , we make 
use of the following result: 



n=0 v ' 



(l-p) N - n R L , n x 



n(n-l) 2n(N-n) (N - n)(N - n + 1) 



p(l -p) 



(1-P) S 



(32) 



Rl, 



Eft 



(28) 



Once the Rl,u is determined, Rl(p), the correspond- 
ing function in the "canonical" percolation ensemble, 
is calculated from Eq. (|T5|), with the binomial coeffi- 
cients for large N being calculated by iterative multi- 
plication uTj. The estimates prg, P0.5, and p c - c for the 



The above three results, along with Eq. (|l5|), demon- 
strate further the advantage of calculating Rl (p) through 
the microcanonical Rl.u- quantities like Rl{p) and 
R^ (p) can be calculated exactly at all p, while p av can be 
found without introducing any error through numerical 
integration. 

In the more familiar binary search method for find- 
ing Rl(p) , p is increased or decreased to narrow the 
bounds on one's estimate of position of the percolation 
point for a given realization of the disorder. This search 
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8 
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n fi1 37fi5fi 


58148fifi 


5842394 


60831 4< 2 ' 


5899755 


101 1925 


532998 


16 


20.0 x 10 9 


0.6069022 


0.5887819 


0.5898858 


0.598828 (2) 


0.5920104 


0.0633761 


0.518117 


32 


20.0 x 10 9 


0.6016319 


0.5914246 


0.5918352 


0.594825 (2) 


0.5926026 


0.0387203 


0.509535 


64 


30.0 x 10 9 


0.5981485 


0.5923179 


0.5924657 


0.593413 (2) 


0.5927391 


0.0233379 


0.504890 


128 


10.0 x 10 9 


0.5959837 


0.5926087 


0.5926613 


0.592952 (2) 


0.5927577 


0.0139703 


0.502476 


256 


40.0 x 10 s 


0.5946742 


0.5927013 


0.5927208 


0.592808 (2) 


0.5927536 


0.0083343 


0.50124 



TABLE I: Various estimates of p c from the simulations (L > 7), where Ns is the number of samples, and from exact expressions 
(L < 7). The cell-to-cell estimate is Po-c f° r superscript (1) and p'-c f° r superscript (2). Errors in the numerical results are 
generally in the last digit quoted. 



process is stopped after some number m of iterations, 
typically about 15, giving to a resolution of 2~ m on the 
estimate for p c . Thus, Rl(p) is evaluated at only a fi- 
nite set of points, and this adds some uncertainty to 
the calculation, beyond the basic statistical error. Given 
that our microcanonical method is also much faster than 
binary search due to its efficient cluster-merging and 
percolation-checking, there seems no reason to use other 
methods. 



IV. RESULTS OF SIMULATIONS 

Simulations were carried out for L = 7, 8, 16, 32, 64, 
128, 256. The results are given in Table |, along with 
exact results from exhaustive enumeration of states for 
small systems with L < 7. The polynomials from which 
the exact results are derived are listed in the appendix. 
We also conducted some simulations at L = 7 to com- 
pare numerical and exact results, and the agreement was 
found to be perfect within the statistical accuracy of the 
simulations. The pseudo-random number generator used 
for the simulations was the four-tap feedback generator 
known as R9689 or gf sr4 |||. 

Error analysis for the simulation data indicates that 
the estimates of p c are accurate to about six figures, and 
the values of R(p c ) are accurate to four or five figures, as 
indicated the table. We also simulated 2.4 x 10 7 samples 
for a system of size L = 512, but the statistical accu- 
racy of the results was insufficient to add anything to the 
present analysis. 

Consider the results for the estimate p max , whose con- 
vergence is non-monotonic. Its value starts below p c for 
small systems, then goes above p c as lattice size passes 
though L ~ 100, and presumably converges to p c from 
above as L — > oo. Indeed, according to Eq. (^4), p max 
has two correction terms with closely-spaced scaling ex- 
ponents — lo — 1/V ~ —1.65 (using the value of lo from 



below) and —1 — 1/V = —1.75; it appears that these 
terms contribute more or less equally in the range of sys- 
tem sizes that we are considering. The observed behavior 
is consistent with Eq. ( p4| ) if 62 < 0, C2 > 0, and the two 
are are roughly comparable in magnitude. (Note that 
03 < because /o(x) is at a maximum at x = 0). It is 
not possible to fit the exponents of Eq. (|24| ) reliably to 
these data. 

The rest of the estimates are all monotonic, and lead 
to reasonably straight lines when viewed on a logarithmic 
plot of \p C st — Pc\ vs. L, reflecting the leading power-law 
behavior. To provide a more sensitive representation of 
our data, we calculate successive slopes between pairs 
of points for systems of size and L, giving values 

- [log |pest (L) -p c | - log|pe S t(^i) -p c |]/log2 for the 
leading exponent for the various estimates p est (L). For 
these calculations we used the value p c — 0.5927462 given 
in Ref. [l^, which is consistent with the data presented 
here, but of somewhat higher precision than these data 
would yield. 

In Fig. Fy we show the plots of the successive slopes 
(2) 

for estimates po.5(L), p v c _ c {L), and p av (L), as a function 
of L~ u with lo = 0.9. According to Eqs. (§|) and @, 
both of these estimates should converge with a leading 
exponent of —1 — 1/V and a next-order term of order 
— 1 — l/V — ui, which implies that the successive slopes 
should fall on a straight line when plotted as a function 
of L~ u , with intercept of 1.75. This behavior is indeed 
seen in Fig. [l], with measured intercepts of 1.754 and 
1.763 respectively. (Note that agreement is not highly 
sensitive to the value of lo ; if the data were plotted as a 
function of L~ l , the fit to linearity would not be much 
worse.) 

The successive slopes for p av (L) do not fall on as good a 
straight line as the other estimates, presumably because 
the exponents — lj — 1/v and — 1 — 1/V of the two leading 
terms in the convergence are closely spaced (see Eq. (^0|) ) . 
Extrapolation to L = 00 is still possible however, and we 
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FIG. 1: Values of the pairwise slopes — [log |p cst (L) — p c \ — 
log |p cst (|L) — p c |] /log 2, plotted as a function , for the 

estimates po.s(L) (circles), pav(L) (triangles), and P^1 C {L) 
(squares) with p c = 0.5927462. 



find an intercept at 1/v+u) = 1.65±0.02, clearly different 
from the value of 1.75 for the other estimates, implying 
uj = 0.90 ± 0.02, the figure we have used above. The 
error bars are smaller than the size of the symbols for 
all L except L — 256, as shown in the plot. Note that 
although uj is used in the abscissa of Fig. [l], its value has 
little effect on the determination of u> from the intercept 
of p a v(i)- 

We can also compare the coefficients for the leading 
behavior of estimates to their predicted values. For ex- 
ample, the predicted value of the leading coefficient for 
Po.5, Eq. (§§), is bo/a! = 0.423, using values of b = 0.322 
(see below) and a\ — 0.765 []|, [l2|]. This compares favor- 
ably with the value measured here of 0.436. And for 
pfl c , the coefficient from @, VM 1 - 2 _1/l/ )], should 
have a value of 1.04, which compares favorably with the 
measured value 1.02. Note that if we take the linear com- 
bination of these two estimates (whose finite-size correc- 
tions are opposite in sign), (po.5 + cxp c _ c )/(l + a) where 
a = 1 — 2~ 1 /' y , the leading correction terms are pre- 
dicted to cancel one another and the combination should 
have leading scaling of = L~ 2 - 65 . And in- 

deed this combination is seen to converge very quickly 
in our numerical results, with values 0.592698, 0.592739, 
0.592745, and 0.592746 for L = 32, 64, 128, and 256 re- 
spectively. The plot of these figures vs. L~ 265 given in 
Fig. H shows linear behavior as expected, with an inter- 
cept at p = 0.5927464(5), consistent with the best current 
figure of p c = 0.5927462(1) ju|. Thus, by taking a combi- 
nation of estimates, we can improve the convergence rate 
for the open system to the point where it becomes com- 
petitive with that of the periodic system, in which the es- 
timate with the best convergence has exponent — ^ [jl6| . 
This cancellation of leading order corrections between the 
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L~ 2 ' 65 

FIG. 2: Plot of (po.s + apfl c )/(l + a) vs. L" 2 ' 65 where a = 
1 - 2~ 1/ ". 



two terms is expected to be universal. 

The results for prg and the standard deviation a = 
\J ((Ap) 2 ) = (p 2 ) — (p) 2 converge to the infinite sys- 
tem value with exponent \jv = —0.75 as shown in Fig. 
H where we plot the successive slopes as a function of 
1/L. (The convergence of the standard deviation of the 
distribution is predicted in Ref. [l2| and also follows from 
the equations of Section [nj). The intercepts are at 0.749 
and 0.763 respectively. We also show in Fig. || a plot of 
the type introduced by Stauffer, in which prg is shown 
as a function of a, allowing extrapolation to infinite sys- 
tem size to be carried out without knowledge of the value 
of v jl]. This plot shows nearly linear behavior, with the 
last three points (L = 64, 128, and 256) well fit by the line 
Prg = 0.5927465 + 0.231512cr with R 2 = 0.9999979. The 
intercept is in excellent agreement with the known value 
of p ci although this agreement is perhaps fortuitous, con- 
sidering the slow convergence of the RG estimate. 

The last column of Table || gives the crossing prob- 
ability at p c . The considerations in Section O imply 
that Rl(Pc) ~ 1/2 + bo/L with no contributions from 
the irrelevant scaling variable || |l2|, |2(| , and indeed 
an analysis of this data shows good agreement with the 
behavior R L (p c ) = 0.320/L - 0.44/L 2 + yielding 
&o = 0.320 ± 0.001. This is nearly identical to the value 
0.319 given in Ref. ^| (where larger systems, but with 
lower statistics, were generated) and the value 0.31 ±0.01 
of Ref. [I§ 



V. CONCLUSIONS 

We have studied the finite-size scaling of estimates of 
the percolation threshold p c derived from the crossing 
probability Rl{p) for site percolation on the square lat- 
tice. Our numerical results confirm that different esti- 
mates converge to p c with a variety of scaling exponents 
as predicted by the scaling theory developed in Refs. [jj, 
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FIG. 3: Values of the pairwise slopes of —(log |prg(£) — Pc \ — 
log|p RG (L/2) -Pc|)/ log 2, with p c = 0.5927462 (triangles), 
and of the standard deviation of the distribution of p (circles), 
plotted as a function L~ . The lines are fit through all the 
points. 



|10| , and |12|. In particular, we have shown that the average 
threshold estimate p av converges with a nontrivial expo- 
nent i J ~ u ~ 1 l lJ whose origins lie in the irrelevant variables 
in the renormalization group treatment of the problem, 
and our results for this estimate provide us with a direct 
measurement of that exponent. We find w = 0.90 ± 02, 
somewhat higher than the value of about 0.85 found pre- 
viously Q but consistent with the (wide) bounds set by 
series studies Jl4|, Our result is in good agreement 
with a renormalization-group result of u> = 0.915 found 
by Burkhardt |27], ||^] . More extensive simulations could 
be done to give to to higher precision. 

Of the estimates considered here, p av is the only one 
that shows the effect of the irrelevant exponent clearly 
The estimates po.5 and p c -c are confirmed to converge as 
L~ 1 ~ 1 /", as proposed previously ||. The maximum es- 
timate p m ax is found to exhibit non-monotonic behavior, 
which can be explained by competition between correc- 
tion terms with closely spaced exponents —1 — \ jv and 

— LO — 1/v. 

The numerical results reported here were found using 
a microcanonical simulation method which allows one to 
calculate Rl(p) easily for any p |l7]]. The various 
estimates can then be found quickly to any desired degree 
of precision by applying appropriate formulas, Eqs. (^]- 
0). This method proves particularly advantageous for 
the estimate p av , since this estimate depends on knowing 
Rl (p) for all values of p, the determination of which by 
most other methods requires a great deal of work. From 
the microcanonical data, p av can be found without any 
calculational bias using Eq. (|30|). 

Having characterized the convergence rates of our var- 
ious threshold estimates, one can go back to older liter- 
ature and find many instances where an anomalous rate 
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FIG. 4: "Stauffer plot" of p RG vs. a = ^/((Ap) 2 )- the line is 
fit through the leftmost three points and its equation is given 
in the text. 



of convergence was seen but not fully understood or rec- 
ognized. For example, Reynolds et al. g] derived an esti- 
mate of p c which is equivalent to our p av from numerical 
data for the one-way crossing probability which they de- 
noted Ri. Assuming this estimate to scale as Eq. (||), 
they plotted their results against L~ x l v (their Fig. 13). 
The resulting plot is seen to fall on a nearly vertical line, 
consistent with higher-order behavior. Fitting their data 
with the supposed L~ x l v scaling, they deduced a best 
estimate of p c — 0.5931 in the large system-size limit. If 
however one assumes instead the L^^^ 1 / u scaling derived 
here, the intercept of their data becomes p c ~ 0.5927, 
which is much closer to the current best estimate of this 
quantity. 

Yonezawa et al. H plotted a quantity essentially equiv- 
alent to our po.5 as a function of L~ x / U , and found 
apparent agreement with this assumed L-dependence 
(their Figs. 7 and 8). The expected behavior 
is evidently too weak to be distinguished from L~ x l v 
within the errors on their numerical data. Similarly, 
Hu et al. believed the cell-to-cell estimate p c 1 } c to 
be insensitive to L\ again, the precision of their work 
did not allow them to observe the higher-order scaling 
predicted by Eq. @. 

It should be emphasized that the convergence behav- 
ior discussed here is specific to a two-dimensional system 
with a square or rhomboidal open boundary, with cross- 
ing defined as a path from one specified side to its oppo- 
site. Because of the symmetry of this system, some terms 
cancel out, allowing various higher-order corrections to 
become dominant. For different boundaries (such as rect- 
angular ones), and for higher dimensional systems, the 
behavior will generally be different. In those cases, most 
of our estimates will scale as the conventional L -1 /", ex- 
cept perhaps the estimates p c - c andp^ c . (To employ the 
latter, R c must be known, but we have exact values only 
for rectangular and conformally related two-dimensional 
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Rl(j>) 

3p 3 + 4p 4 - 6p 5 - 9p 6 + Up 7 - 6p 8 + p 9 

4p 4 + 12p 5 - 6p 6 - 28p 7 - 22p 8 + 48p 9 + 66p 10 - 108p n + 10p 12 + 44p 13 - 20p 14 + p 16 

5p 5 + 24p 6 + 12p 7 - 62p 8 - 92p 9 - 41p 10 + 274p 41 + 42p 12 + 474p 13 - 1336p 14 + 172p 15 - 197p 16 + 4791p 17 - 9015P 1 

+ 8013p 19 - 4261p 20 + 1559p 21 - 450p 22 + 103p 23 - 15p 24 + p 25 
6p 6 + 40p 7 + 60p 8 - 80p 9 - 248p 10 - 276p n + 201p 12 + 944p 13 - 298p 14 + 2392p 15 - 2420p 1B + 548p 17 - 24848p 18 

+ 38688p 19 - 7540p 20 + 5 1564p 21 - 117312p 22 - 1333 12p 23 + 588639p 24 - 6084 64p 25 + 68 3 62p 26 + 4203 96p 27 

- 455910p 28 + 235816p 29 - 62454p 30 + 3200p 31 + 3212p 32 - 1024p 33 + 120p 34 - p 36 

7p 7 + 60p 8 + 150p 9 - 18p 10 - 490p n - 885p 12 - 318p 13 + 1464p 14 + 3056p 15 - 1586p 16 + 5584p 17 - 6520p 18 + 43150p 

- 153589p 20 + 128 5 04p 21 - 407257p 22 + 1278288p 23 - 1243 193p 24 + 2 195 3 74p 25 - 69 8 3 6 30p 26 + 8271536p 27 

- 6990669p 28 + 17741331p 29 - 11344431p 30 - 55 2 9 4 9 29p 31 + 91642905p 32 + 67152194p 33 - 3742 5 5 5 72p 34 
+ 557174473p 35 - 463108229p 36 + 22 5 3 3 8 9 48p 37 - 47135360p 38 - 129 5 691p 39 + 11168 8 48p 40 - 1724804p 41 

- 1067305p 42 + 68 93 18p 43 - 196565p 44 + 34 8 48p 45 - 4391p 46 + 422p 47 - 28p 48 + p 49 

TABLE II: Exact results for Rl{p) expressed as polynomials in p, for L = 2 . . . 7. 



systems p9|.) Study of these other systems is a subject 
for future research. 

Another approach to measuring p c is to use periodic 
rather than open boundary conditions. A partially peri- 
odic system in two dimensions is a cylinder, and crossing 
in this system was studied in Ref. [l2]. The fully peri- 
odic rectangle is a torus, and the criterion of crossing 
is replaced by criteria involving the different topologi- 
cally distinct ways in which clusters can wrap around 
the boundaries [j3(J. (Some authors (ll], |24|, |3lJ have also 
considered the percolation criterion in which a cluster has 
the full dimension of the lattice along at least one axis 
but does not necessarily wrap around.) In Ref. [l7| we 
showed that many estimates of p c on the torus converge 
a factor L faster than estimates for the open square — 
some converging as fast as L^ 2 ^ x l v . 

In conclusion, it is clear that the convergence of esti- 
mates for the critical occupation probability p c in per- 
colation systems is highly dependent upon the nature of 
the estimate, as well as the shape and boundary condi- 
tions of the system, and that the shrewd use of this fact 
can allow one to make very accurate estimates of p c and 
scaling exponents. 



APPENDIX A: 



EXACT ENUMERATION 
RESULTS 



In Table |j] we give the exact expressions for the cross- 
ing probability function Rl{p) for site percolation on a 
square lattice of size L x L, for crossing from one given 
side of the square to the opposite side (such as left to 
right). The results for L = 2 to 5 were given previously 
by Reynolds et al. Those for L — 6 and L = 7 were 
calculated previously for the work reported in Ref. ||, 



but reported in a different format (results were given for 
R' L {p) rather than Rl(p) itself). From the results here, 
one can with reasonable ease calculate the various esti- 
mates of p c given in Table || using a symbolic manipu- 
lation program such as Maple or Mathematica. A file 
containing these polynomials in forms readable by such 
programs is available by email from the authors. 

An alternative way to represent these results is as a se- 
ries in p n q N ~ n where q = 1 — p and N = L 2 . The trans- 
formation can be achieved by substituting p — » 1/(1 + r), 
multiplying by (1 + r) N p N , expanding, and replacing 
r ~~ * q/p again. The results are given in Table Jffj]. 

This is also the form that Reynolds et al. [g| used in 
their series for i?2 through R§. From the present point 
of view, these series are interesting because they are in 
precisely the form of Eq. (|l^), so that the coefficient cl,« 
of p n q N ~ n above is related to the microcanonical crossing 
probability R n simply by 



Rl,u — 



CL,r, 
(N\ 



c L , n n\{N - n)\ 
N\ 



(Al) 



That is, CL^ n represents the number of configurations 
with n occupied sites that satisfy the crossing criterion, 
out of a total of (^) possible configurations of the n oc- 
cupied sites among the N = L 2 sites of the lattice. For 
example, of the 16!/(6!10!) = 8008 possible configura- 
tions of 6 occupied sites on a 4 x 4 lattice, exactly 390 are 
percolating by crossing in one direction (from the third 
term in Ri(p,q)), yielding a microcanonical probability 
#4,6 = 390/8008 = 0.048701 . . . Likewise, for n = 16 
occupied sites on the 4x4 system, there is exactly one 
percolating system out of one total system. 

Thus, the polynomials given in Table III represent the 
microcanonical Rl. u for L up to 7. 
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